The GNSS navigation problem is the problem of estimating the position of a GNSS user by means of the information provided by the GNSS signal as received by a GNSS user receiver.
There are several standard GNSS navigation techniques, the most common of which is absolute navigation. In absolute navigation, the navigation system computes its absolute position and velocity with no more information than that contained in the GNSS satellite signals tracked by the user's receiver, by means of so-called pseudo-range and Doppler measurements. For that purpose, in order to be able to obtain its position some additional parameters must be estimated, as the receiver clock bias and drift. If more than one satellite constellation is processed (e.g. GPS+GLONASS), it is also necessary to add an inter-system bias that accounts for the different synchronization of the satellite clocks. Position, velocity and the auxiliary parameters are simultaneously estimated in a common Kalman filter.
GNSS navigation can also profit from the use of the data from an inertial measurement unit, or IMU, which usually consists of a 3-axis accelerometer and a 3-axis gyroscope. The GNSS+IMU tightly coupled Kalman filter substitutes the propagation model of the GNSS-only Kalman filter by a mechanization stage, which uses the data from the IMU to predict the filter state throughout the time period between successive GNSS data epochs. This hybrid filter requires the inclusion of new parameters in the Kalman state vector, namely, the accelerometer and gyro biases, plus the attitude matrix variables.
Finally, it is possible to perform absolute navigation with centimeter-level errors by using a technique called PPP (Precise Point Positioning). PPP processes pseudorange and carrier-phase measurements from a single user receiver, using detailed physical models and corrections, and requires precise satellite orbits and clocks products (e.g. through real-time corrections to the navigation message). PPP is different from other precise-positioning approaches like RTK in that no reference stations are needed in the vicinity of the user receiver. The only observation data that must be processed are measurements from the user receiver. In terms of the Kalman filter, the main difference of PPP with respect to the standard GNSS navigation is that a carrier-phase ambiguity parameter must be estimated for each satellite.
The present invention is applicable in the GNSS standalone, GNSS+IMU and PPP cases.
For many applications of GNSS navigation a certain level of accuracy is not enough. In addition, the system must also provide reliability. For instance, in a tolling system users should always be charged for all and only the toll roads they use (no more, no less). In other fields of application, as aviation, the consequences of getting a wrong position beyond certain limits and not being aware of it are much more dramatic. This can be solved, in the case of a GNSS-based system, by adding a method of navigational integrity. Integrity provides the necessary guarantee that users are where they are supposed to be (i.e. where the navigation system says they are). A system can preserve integrity and at the same time be inaccurate, as long as the system itself is aware of the lack of accuracy. In such a case, the problem is that the system may be unavailable more often than desired. A perfect system (from the point of view of its reliability) would exhibit 100% availability and 100% integrity. Unfortunately, perfect systems do not exist, and in practice there is always a trade-off between both goals.
Many methods of integrity use the concept of Protection Level (PL), which intends to provide bounds to the GNSS positioning errors up to certain level of confidence. For instance, for a given confidence level 1−α, a valid Protection Level is a positive number B such that:P(|δ|>B)<α
where P is the probability operator and |δ| is the norm of the position error. Usual values of the confidence level α range from 10−4 to 10−7, and depend on the requirements of the GNSS application. For instance, α=10−7 appears frequently in the aviation context. Note that there are many possible values that fulfil the Protection Level condition, but it is important that the method provides fairly small ones, so that the availability is not degraded. In the case of GNSS navigation, the concept of error bound is usually particularized to a subset of the position coordinates, especially the vertical component δU (e.g. civil aviation) or the horizontal ones, in the North and East directions (δE,δN), (e.g. electronic road fee collection).
This concept of Protection Level arose as the core of the GNSS Integrity concept that was developed for Satellite Based Augmentation Systems (SBAS), such as the American WAAS or the European EGNOS among others, and has been applied specifically to those systems as well as to Ground Based Augmentation Systems (GBAS) such as LAAS. Both SBAS and GBAS are systems that provide real time corrections and integrity information applicable to GNSS observations. Civil aircrafts equipped with special GNSS navigation units can receive those corrections and apply them to their GNSS observations, as well as the integrity information, in the form of statistical bounds to the observation errors that remain after applying the corrections. Thus the on-board GNSS unit can achieve a more accurate estimate of its position (thanks to corrections) and, moreover, can compute a Protection Level for the remaining position error.
Besides, autonomous methods for computing a Protection Level have been defined in the frame of the so-called RAIM methods (Receiver Autonomous Integrity Monitoring). By autonomous it is meant that they do not depend on corrections or any extra information coming from an augmentation system such as a SBAS or a GBAS. The RAIM concept aims to provide an integrity layer to the GNSS navigation process, implementing techniques for detecting and isolating faulty measurements (that is, measurements with excessive error) along with the mentioned Protection Levels for statistically bounding the position estimation error. Such PL computation methods are, however, difficult to justify from a theoretical point of view, since they rely on hypotheses that rarely hold in real world.
Patent document EP2113786 discloses an autonomous method to provide Protection Levels for GNSS positioning based on navigation residuals, which does not make any hypothesis about the statistical behaviour of the errors of the individual measurements, but it is based on only one hypothesis (isotropy) that such measurement errors combine in an error vector which can point in any direction (of the measurement space) with the same probability.
All the methods to compute Protection Levels that have been mentioned above only apply to the least-squares estimation. A method to compute Protection Levels for a Kalman solution faces new requirements:                The filter makes use of different types of measurements (e.g. pseudoranges, Doppler and carrier-phase measurements in the GNSS case).        The filter combines observations from different epochs, giving more weight to the most recent ones. In addition, the measurement noise is not white noise, as the Kalman filter usually assumes. Instead, measurements performed in successive epochs show a temporal correlation (coloured noise) that may have an important (negative) impact in the accuracy of the Kalman solution. So, it is crucial that the method to compute the Protection Levels takes into account the effects of temporal correlations.        The Kalman solution has errors coming from the propagation of the solution between different updates, which should also be quantified.        The size of the Kalman state vector may become quite big, but usually only a few of the estimated parameters are interesting.        
Multivariate t-distributions can be used to model different error distributions. The family is parameterized by a scalar Nm, which is called the number of degrees of freedom, and a covariance matrix Rm, which gives an estimation of the size of the errors. Thus, the distribution of an error vector θ of dimension d reads:
                              θ          ~                                    𝒯              d                        ⁡                          (                                                N                  m                                ,                                  R                  m                                            )                                      =                                            Γ              ⁡                              [                                                      (                                          n                      +                      d                                        )                                    /                  2                                ]                                                                    Γ                ⁡                                  (                                      n                    /                    2                                    )                                            ⁢                              N                m                                  d                  /                  2                                            ⁢                              π                                  d                  /                  2                                            ⁢                                                                                      R                    m                                                                                      1                  /                  2                                                              ·                                    [                              1                +                                                      1                                          N                      m                                                        ⁢                                      θ                    T                                    ⁢                                      R                    m                                          -                      1                                                        ⁢                  θ                                            ]                                                      -                                  (                                                            N                      m                                        +                    d                                    )                                            /              2                                                          [                  Eq          .                                          ⁢          1                ]            
The t-distribution is heavy tailed and it provides a more realistic representation of the measurement and navigation errors than the normal distribution. It has been used in order to improve the robustness of navigation filters against outliers. It provides a good model for the measurement errors, which is used in a Bayesian framework (see, for instance, “A New Approach for Improving Reliability of Personal Navigation Devices under Harsh GNSS Signal Conditions”, Anup Dhital, Jared B. Bancroft and Gérard Lachapelle, Sensors, November 2013, 13, 15221-15241).
The following paragraphs provide a summary of the Kalman filter and sets up the notation that will be used later. The Kalman filter estimates the value at successive epochs of a set of parameters (state vector), by using the information of the measurements that arrive sequentially. The measurements are processed in a sequential way, providing a new solution every time a new set of observations are received. The filter proceeds by updates of its main variables:                xk, the state vector, or vector of parameters defining the state at time or epoch tk;        Pk, the covariance matrix of the estimated parameters at time tk;        
It performs two different types of updates:
A. Time Update (“Propagation”):
The filter propagates the state vector and covariance matrix from epoch k−1 to epoch k; their value after this update is usually noted xk−, Pk−. The filter models the stochastic behaviour of the parameters between successive epochs. In the linear setting—which is enough for most purposes—the state and covariance propagation is given by:xk−=Fkxk−1 Pk−=FkPk−1FkT+Qk where:                Fk is the transition matrix from time tk−1 to tk; and,        Qk is the covariance matrix of the random process (“process noise”) that drives the evolution of the parameters of the state vector.        
In the hybrid filter, the propagation of the state and its covariance matrix is performed by a mechanization algorithm that processes the IMU data received. For most purposes, this update can be considered as equivalent to the linear update.
B. Measurements Update:
The second update uses the information provided by the measurements to correct the state vector and covariance matrix, and gives their final value at epoch tk, noted xk and Pk, respectively. The value of the measurements vector at epoch k can be written aszk=h(xk)+vk 
Function h gives the dependence of the measurements on the parameters, and vk is the measurement noise, assumed to have a covariance matrix Rk. The estimation is based on the linearized expression:zk=HkTxk+vk 
The matrix Hk is usually known as the design matrix (also called observation matrix). The measurements update can be performed through the following steps:                The measurements are reconstructed to obtain the residuals vector yk;        The new covariance matrix Pk is computed as follows:Pk=(HkTWkHk+(Pk−)−1)−1         where Wk=Rk−1 is the weight matrix of the measurements at epoch k.        The normal equations are solved:Δxk=PkHkTWk·yk         so that xk=xk−+Δxk.        
There are other mechanizations of the sequential filter, which are mathematically equivalent, but the one described above is preferred to present the integrity algorithm. For reasons of clarity, in the rest of the document the ubiquitous subscript k will be dropped. The design and weight matrix will be decomposed in the blocks corresponding to each measurement type m, noted as Hm, Wm respectively, as shown in the following diagrams:
      H    =          (                                                  H                              m                ⁢                                                                  ⁢                1                                                                          ⋮                                                              H              mk                                          )        ,      W    =          (                                                  W                              m                ⁢                                                                  ⁢                1                                                          …                                0                                                ⋮                                ⋱                                ⋮                                                0                                …                                              W              mk                                          )      
Finally, it is important to discuss the influence of the temporal correlations of the measurements in the performance of the filter. The Kalman filter as presented above is optimal under the assumption that the measurements taken at different epochs are uncorrelated. However, this is seldom the case. Usually, the measurement process is repeated at each epoch and give a constant number of measurements (say w), with an associated noise vector vk of size w. The noise is assumed to have zero mean: E(vk)=0. The temporal correlation matrices can be defined as B(k,l)=E(vkvlT). The noise process is called stationary if B(k,l)=B(k+s,l+s)=B(k−l). In the case of the white noise, the measurements at different epochs are not correlated, and then B(k,l)=0 if k≠1. In most cases the correlation matrix B is diagonal, and each component of the vector vk can be analysed as a one-dimensional process. In one dimension we can define the autocorrelation function:
      R    ⁡          (              k        ,        l            )        =            B      ⁡              (                  k          ,          l                )                                      B          ⁡                      (                          k              ,              k                        )                                    1          /          2                    ⁢                        B          ⁡                      (                          l              ,              l                        )                                    1          /          2                    
If the process is stationary, the autocorrelation function is simpler and may be written as:R(k,l)=ρ(k−l)
The function ρ verifies the following inequalities:0≤ρ(h)≤1
Then, the diagonal correlation matrix B of a stationary process will be:
            B      ⁡              (        h        )              =          (                                                                  ρ                1                            ⁡                              (                h                )                                                          …                                0                                                ⋮                                ⋱                                ⋮                                                0                                …                                                              ρ                w                            ⁡                              (                h                )                                                        )        ,      0    ≤                  ρ        i            ⁡              (        h        )              ≤    1  
In the simplest case, a Gauss-Markov process, the autocorrelation function is given by:ρ(h)=φ|h|, 0≤φ≤1
For values of ρ(h) significantly above 0, the filter performance is clearly degraded. So, it is very important to take into account the effect of correlations to compute bounds for the Kalman solution errors. In the method of the present invention, when for a measurement type m the temporal correlations are considered, we will denote by ρm the associated correlation matrix B(1) of the measurement noise. This matrix will be given as an input to the method. Its values will be derived either from the knowledge of the measurement device characteristics or by any independent method that estimates them from the evolution of the measurements and their residuals.
As far as the present invention is concerned, the actual form of the non-linear GNSS observation equation or how the design matrix H is derived from it, are not relevant topics, so the present document does not go deeper into such details (which, on the other hand, are of standard use in GNSS Kalman filter navigation and can be learned from many GNSS literature sources as, for instance, “Global Positioning System: Theory & Applications”, Bradford W. Parkinson & James J. Spilker (editors), 1996).